#ifndef YPL_DFT_H
#define YPL_DFT_H

#include "Array1D.h"
#include "Array2D.h"
#include "SComplex.h"

namespace ypl
{

class DFT
{
public:
    /*! 傅立叶变换的方向 */
    enum Direct
    {
        Forward,  ///< 正变换
        Backward  ///< 逆变换
    };


    /*************************
     * 采用复数数组的接口
     *************************/

    /*!
     * @brief 一维离散傅里叶变换（复数输入）
     */
    void operator()(const Array1D<Complex> &aIn,
                    Array1D<Complex> &aOut,
                    Direct direct = Forward);
    /*!
     * @brief 一维离散傅里叶变换（实数输入）
     */
    void operator()(const Array1D<Real> &aIn,
                    Array1D<Complex> &aOut,
                    Direct direct = Forward);
    /*!
     * @brief 二维离散傅里叶变换（复数输入）
     */
    void operator()(Array2D<Complex> &aIn,
                    Array2D<Complex> &aOut,
                    Direct direct = Forward);
    /*!
     * @brief 二维离散傅里叶变换（实数输入）
     */
    void operator()(const Array2D<Real> &aIn,
                    Array2D<Complex> &aOut,
                    Direct direct = Forward);

    /*****************************
     * 采用简单复数数组的接口
     *****************************/

    /*!
     * @brief 一维离散傅里叶变换（复数输入）
     */
    void operator()(const Array1D<SComplex> &aIn,
                    Array1D<SComplex> &aOut,
                    Direct direct = Forward);
    /*!
     * @brief 一维离散傅里叶变换（实数输入）
     */
    void operator()(const Array1D<Real> &aIn,
                    Array1D<SComplex> &aOut,
                    Direct direct = Forward);
    /*!
     * @brief 二维离散傅里叶变换（复数输入）
     */
    void operator()(Array2D<SComplex> &aIn,
                    Array2D<SComplex> &aOut,
                    Direct direct = Forward);
    /*!
     * @brief 二维离散傅里叶变换（实数输入）
     */
    void operator()(const Array2D<Real> &aIn,
                    Array2D<SComplex> &aOut,
                    Direct direct = Forward);

protected:

    void dft(const SComplex *pIn, Integral n,
             SComplex *pOut, Direct direct = Forward);
    void fft(const SComplex *pIn, Integral n,
             SComplex *pOut, Direct direct = Forward);
    void dft(SComplex **pIn, Integral dim1, Integral dim2,
             SComplex **pOut, Direct direct = Forward);
    void fft(SComplex **pIn, Integral dim1, Integral dim2,
             SComplex **pOut, Direct direct = Forward);

private:

    void bitReverse(Integral n, Array1D<Integral> &a);
    void orderCopy(const SComplex *pIn, Integral n,
                   SComplex *pOut,
                   const Array1D<Integral> &aOrder);
    void calcW(Array1D<SComplex> &aW, Integral n,
               Integral m = 0,
               Direct direct = Forward);
    void calcKernel(SComplex *pData, Integral dim,
                    const SComplex *pW);
    void calcKernelDim1(SComplex **pData, Integral dim1, Integral dim2,
                        const SComplex *pW);
    void calcKernelDim2(SComplex **pData, Integral dim1, Integral dim2,
                        const SComplex *pW);
    void divide(SComplex *pData, Integral dim,
                Integral divisor = 0);
    void divide(SComplex **pData, Integral dim1, Integral dim2,
                Integral divisor = 0);
    Integral log2(Integral n);
    bool isPower2(Integral n);
};



/*!
 * @brief 复数序列按指定顺序复制到简单复数序列
 * @param pIn       -in-    原复数序列
 * @param n         -in-    复数序列点数
 * @param psOut     -out-   简单复数序列
 * @param aOrder    -in-    指定新顺序（可选参数），默认为原顺序
 */
inline
void DFT::orderCopy(const SComplex *pIn,
                    Integral n,
                    SComplex *pOut,
                    const Array1D<Integral> &aOrder)
{
    /*! 原序复制 */
    if (aOrder.dim() != n)
    {
        for (Integral i = 0; i < n; ++i)
        {
            pOut[i].real = pIn[i].real;
            pOut[i].imag = pIn[i].imag;
        }
        return;
    }
    /*! 按指定顺序复制 */
    const Integral *pOrder = aOrder;
    for (Integral i = 0; i < n; ++i)
    {
        Integral iNew =  pOrder[i];
        pOut[iNew].real = pIn[i].real;
        pOut[iNew].imag = pIn[i].imag;
    }
}




}  // namespace ypl

#endif
